ON THE NUMERICAL EVALUATION OF FREDHOLM DETERMINANTS 

FOLKMAR BORNEMANN* 

Abstract. Some significant quantities in mathematics and physics are most naturally expressed as the 
Fredholm determinant of an integral operator, most notably many of the distribution functions in random 
matrix theory. Though their numerical values are of interest, there is no systematic numerical treatment 
of Fredholm determinants to be found in the literature. Instead, the few numerical evaluations that are 
available rely on eigenfunction expansions of the operator, if expressible in terms of special functions, 
or on alternative, numerically more straightforwardly accessible analytic expressions, e.g., in terms of 
Fainleve transcendents, that have masterfully been derived in some cases. In this paper we close the gap 
in the literature by studying projection methods and, above all, a simple, easily implementable, general 
method for the numerical evaluation of Fredholm determinants that is derived from the classical Nystrom 
method for the solution of Fredholm equations of the second kind. Using Gauss-Legendre or Clenshaw- 
C5 Curtis as the underlying quadrature rule, we prove that the approximation error essentially behaves 

^ like the quadrature error for the sections of the kernel. In particular, we get exponential convergence 

^~5 for analytic kernels, which are typical in random matrix theory. The application of the method to the 

_^ distribution functions of the Gaussian unitary ensemble (GUE), in the bulk and the edge scaling limit, 

-yi is discussed in detail. After extending the method to systems of integral operators, we evaluate the two- 

point correlation functions of the more recently studied Airy and Airyj processes. 



00 

o 

O 
(N 



Key words. Fredholm determinant, Nystrom's method, projection method, trace class operators, ran- 
dom matrix theory, Tracy-Widom distribution. Airy and Airy^ processes 



(-H AMS subject classifications. 65R20, 65F40, 47G10, 15A52 

a 

r^ 1. Introduction. Fredholm's (1903) landmark paper on linear integral equa- 

I "1 tions is generally considered to be the forefather of those mathematical concepts that 
finally led to modern functional analysis and operator theory — see the historical ac- 

\^ counts in Kline (1972, pp. 1058-1075) and Dieudonne (1981, Chap. V). Fredholm was 

J^ interested in the solvability of what is now called a Fredholm equation of the second 

■^ kind, 

u{x)+z K{x,y)u{y)dy = f{x) {x e {a,b)), (1.1) 

o 

QQ and explicit formulas for the solution thereof, for a right hand side / and a kernel K, 

f^ both assumed to be continuous functions. He introduced his now famous determi- 

j^ nant 

X 

c^ k=o 



H ^^^^ ^J^.^.Ja '"la '^''^ ^^^^'" ^''^^M=l "^^^ ■ ■ ■ '^^"' ^"'^^ 



which is an entire function of z G C, and succeeded in showing that the integral 
equation is uniquely solvable if and only if d{z) 7^ 0. 

Realizing the tremendous potential of Fredholm's theory Hilbert started work- 
ing on integral equations in a flurry and, in a series of six papers from 1904 to 1910,^ 

*Zentrum Mathematik, Technische Universitat Mtinchen, Boltzmarmstr. 3, 85747 Garching, Germany 
(bornemann@ma.tum.de). Manuscript as of June 24, 2008. 

'Simon (2005, p. VII) writes: "This deep paper is extremely readable and I recommend it to those 
wishing a pleasurable afternoon." An English translation of the paper can be found in Birkhoff (1973, 
pp. 449-465)- 

^Later reproduced as one of the first books on linear integral equations (Hilbert 1912). 
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transformed the determinantal framing to the beginnings of what later, in the hands 
of Schmidt, Carleman, Riesz, and others, would become the theory of compact op- 
erators in Hilbert spaces. Consequently, over the years Fredholm determinants have 
faded from the core of general accounts on integral equations to the historical re- 
marks section — if they are mentioned at all.^ 

So, given this state of affairs, why then study the numerical evaluation of the 
Fredholm determinant d{z)? The reason is, simply enough, that the Fredholm deter- 
minant and the more general notions, generalizing (i.i) and (1.2) to 

u + zAu=f, d{z) = det{I + zA), 

for certain classes of compact operators A on Hilbert spaces, have always remained 
important tools in operator theory and mathematical physics (Gohberg, Goldberg 
and Krupnik 2000, Simon 2005). In turn, they have found many significant appli- 
cations: e.g., in atomic collision theory (Jost and Pais 1951, Moiseiwitsch 1977), in- 
verse scattering (Dyson 1976), in Floquet theory of periodic differential equations 
(Eastham 1973), in the infinite-dimensional method of stationary phase and Feyn- 
man path integrals (Albeverio and Hoegh-Krohn 1977, Rezende 1994), as the tw^o- 
point correlation function of the two-dimensional Ising model (Wilkinson 1978), in 
renormalization in quantum field theory (Simon 2005), as distribution functions in 
random matrix theory (Mehta 2004, Deift 1999, Katz and Sarnak 1999) and com- 
binatorial growth processes (Johansson 2000, Prahofer and Spohn 2002, Sasamoto 
2005, Borodin, Ferrari, Prahofer and Sasamoto 2007). As Lax (2002, p. 260) puts it 
most aptly upon including Fredholm's theory as a chapter of its own in his recent 
textbook on functional analysis: "Since this determinant appears in some modern 
theories, it is time to resurrect it." 

In view of this renewed interest in operator determinants, w^hat numerical meth- 
ods are available for their evaluation? Interestingly, this question has apparently 
never — at least to our knowledge — been systematically addressed in the numerical 
analysis literature.4 Even experts in the applications of Fredholm determinants com- 
monly seem to have been thinking (Spohn 2008) that an evaluation is only possible 



3For example, Hilbert (1912), Webster (1927), and Whittaker and Watson (1927) start with the Fred- 
holm determinant right from the begirming; yet already Courant and Hilbert (1953, pp. 142-147), the 
translation of the German edition from 1931, give it just a short mention ("since we shall not make any 
use of the Fredholm formulas later on"); while Smithies (1958, Chap. V) and Hochstadt (1973, Chap. 7), 
acknowledging the fact that "classical" Fredholm theory yields a number of results that functional ana- 
lytic techniques do not, postpone Fredholm's theory to a later chapter; whereas Baker (1977), Porter and 
Stirling (1990), Prossdorf and Silbermann (1991), Hackbusch (1995), and Kress (1999) ignore the Fredholm 
determinant altogether. Among the newer books on linear integral equations, the monumental four vol- 
ume work of Fenyo and Stolle (1982-1984) is one of the few we know of that give Fredholm determinants 
a balanced treatment. 

^Though we can only speculate about the reasons, there is something like a disapproving attitude 
towards determinants in general that seems to be quite common among people working in "continuous" 
applied mathematics. Here are a few scattered examples: Meyer (2000) writes at the beginning of the 
chapter on determinants (p. 460): "Today matrix and linear algebra are in the main stream of applied 
mathematics, while the role of determinants has been relegated to a minor backwater position." Axler 
(1995) has a paper with the provocative title "Down with Determinants!" and a correspondingly worked 
out textbook on linear algebra (1997). The quintessential book of Golub and Van Loan (1996) on "Matrix 
Computations" does not explicitly address the computation of determinants at all, it is only implicitly 
stated as part of Theorem 3.2.1. Higham (2002) writes at the beginning of Section 14.6: "Like the matrix 
inverse, the determinant is a quantity that rarely needs to be computed." He then continues with the 
argument, well known to every numerical analyst, that the determinant cannot be used as a measure 
of ill conditioning since it scales as det(aA) = a'" det{A) for a m x wi-matrix A, a G R. Certainly there 
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if either the eigenvalues of the integral operator are, more or less, explicitly known 
or if an alternative analytic expression has been found that is numerically more 
accessible — in each specific case anew, lacking a general procedure. 

The Nystrom-type method advocated in this paper. In contrast, we study a simple 
general numerical method for Fredholm determinants which is exceptionally effi- 
cient for smooth kernels, yielding small absolute errors (i.e., errors that are small 
with respect to the scale det(I) = 1 inherently given by the operator determinant). 
To this end w^e follow the line of thought of Nystrom's (1930) classical quadrature 
method for the numerical solution of the Fredholm equation (1.1). Namely, given a 
quadrature rule 

m .b 

Q{f) = T.^"if{^j)^ / f{^)dx, 
1=1 ■''' 

Nystrom discretized (1.1) as the linear system 

m 

Ui + zY^WjK{xi,Xj)Uj = f{xi) (f = l,...,m), (1.3) 

/=i 

which has to be solved for m, w m(x,) {i — 1, . . . ,m). Nystrom's method is extremely 
simple and, yet, extremely effective for smooth kernels. So much so that Delves and 
Mohamed (1985, p. 245), in a chapter comparing different numerical methods for 
Fredholm equations of the second kind, write: 

Despite the theoretical and practical care lavished on the more com- 
plicated algorithms, the clear winner of this contest has been the 
Nystrom routine with the wz-pornt Gauss-Legendre rule. This rou- 
tine is extremely simple; it includes a call to a routine which pro- 
vides the points and weights for the quadrature rule, about twelve 
Unes of code to set up the Nystrom equations and a call to the rou- 
tine which solves these equations. Such results are enough to make 
a numerical analyst weep. 

By keeping this conceptual and algorithmic simplicity, the method studied in this 
paper approximates the Fredholm determinant d(z) simply by the determinant of 
the m X w2-matrix that is applied to the vector (m,) in the Nystrom equation (1.3): 

dQ{z) = det {6ij + z w;jX(x;,xy))|".^^ . (1.4) 

If the w^eights if, of the quadrature rule are positive (which is always the better 
choice), w^e will use the equivalent symmetric variant 

dQ{z) = det [&i^^zw]I^K{xi,x-;)w]I^Y-^,^^ . (1.5) 



is much truth in all of their theses, and Cramer's rule and the characteristic polynomial, which were 
the most common reasons for a call to the numerical evaluation of determinants (Stewart 1998, p. 176), 
have most righteously been banned from the toolbox of a numerical analyst for reasons of efficiency. 
However, with respect to the infinite dimensional case, the elimination of determinants from the thinking 
of numerical analysts as a subject of computations might have been all too successful. For instance, the 
scaling argument does not apply in the infinite dimensional case: operator determinants det(7 + A) are 
defined for compact perturbations of the identity, which perfectly determines the scaling since, for a 7^ 1, 
B.{1 + A) cannot be written in the form / + A with another compact operator A. (This is because the 
identity operator is noi compact then.) 
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Using Gauss-Legendre or Curtis-Clenshaw quadrature rules, the computational 
cost^ of the method is of order 0{m^). The implementation in Matlab, or Mathe- 
matica, is straightforward and takes just a few lines of code.^ In Matlab: 

function d = DetNystr6m(K,z,a,b,m) 

[w,x] = QuadratureRule(a,b,m) ; 

w = sqrt (w) ; 

[xj,xi] = meshgrid(x,x) ; 

d = det (eye(m)+z* (w' *w) . *K(xi,xj)) ; 

In Mathematica: 

DetNystr6m[K_, z_, a_, b_, m_] : = 

Module {w, x} , {w, x} = QuadratureRule [a, b, m] ; w = Vw ; 
Det [IdentityMatrix[m] + z Outer [Times, w, w] Outer[K, x, x] ] ; 

Strictly speaking we are not the first to suggest this simple method, though. 
In fact, it w^as Hilbert himself (1904, pp. 52-60) in his very first paper on integral 
equations, w^ho motivated^ the expression (1.2) of the Fredholm determinant by es- 
sentially using this method with the rectangular rule for quadrature, proving locally 
uniform convergence; see also Hilbert (1912, pp. 4-10) and, for the motivational 
argument given just heuristically without a proof of convergence, Whittaker and 
Watson (1927, Sect. 11.2), Tricomi (1957, pp. 66-68) (who speaks of a "poetic license" 
to be applied "without too many scruples"). Smithies (1958, pp. 65-68), Hochstadt 
(1973, pp. 243-239), and Fenyo and Stolle (1982-1984, Vol. II, pp. 82-84) — to name 
just a few but influential cases. Quite astonishingly, despite of all its presence as 
a motivational tool in the expositions of the classical theory, w^e have found just 
one example of the use of this method (with Gauss-Legendre quadrature) in an ac- 
tual numerical calculation: a paper by the physicists Reinhardt and Szabo (1970) on 
low-energy elastic scattering of electrons from hydrogen atoms. However, the error 
estimates (Theorem 6.2) that w^e will give in this paper seem to be new at least; we 
will prove that the approximation error essentially behaves like the quadrature error 
for the sections x h^ K{x,y) and y h^ K{x,y) of the kernel. In particular, w^e will get 
exponential convergence rates for analytic kernels. 

Examples. Perhaps the generality and efficiency offered by our direct numerical 
approach to Fredholm determinants, as compared to analytic methods if they are 

'The computational cost of 0{nr^) for the matrix determinant using either Gaussian elimination with 
partial pivoting (Stewart 1998, p. 176) or Hyman's method (Higham 2002, Sect. 14.6) clearly dominates 
the cost of 0(m log m) for the weights and points of Clenshaw-Curtis quadrature using the FFT, as well as 
the cost of 0{m-^) for Gauss-Legendre quadrature using the Golub-Welsh algorithm; for implementation 
details of these quadrature rules see Waldvogel (2006), Trefethen (2008), and the appendix of this paper 
The latter paper carefully compares Clenshaw-Curtis with Gauss-Legendre and concludes (p. 84): "Gauss 
quadrature is a beautiful and powerful idea. Yet the Clenshaw-Curtis formula has essentially the same 
performance for most integrands." 

^The command [w,x] = QuadratureRule(a,b,m) is supposed to supply the weights and points of 
a m-point quadrature rule on the interval [a, fc] as a 1 x m vector w and am x 1-vector x, respectively. For 
Gauss-Legendre and Clenshaw-Curtis, such a Matlab code can be found in the appendix. 

^Fredholm himself does not give the slightest hint of a motivation in his early papers (1900, 1903). He 
apparently conceived his determinant in analogy to similar expressions that his fellow countryman von 
Koch (1892) had obtained for infinite matrices a few years earlier; see Fredholm (1909, p. 95), Dieudonne 
(1981, p. 99), or Pietsch (2007, p. 409). 
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available at all, is best appreciated by an example. The probability £2(0; s) that an 
interval of length s does not contain, in the bulk scaling limit of level spacing 1, 
an eigenvalue of the Gaussian unitary ensemble (GUE) is given by the Fredholm 
determinant of the sine kernel, 

E2(0;s) = det (l - As\i^2^o,s)) . Asw(x) = ^' ^^^^^^^^u{y) dy; 

see Gaudin (1961) and Mehta (2004, Sect. 6.3). Gaudin has further shown that the 
eigenfunctions of this selfadjoint integral operator are exactly given by a particular 
family of special functions, namely the radial prolate spheroidal wave functions with 
certain parameters. Using tables (Stratton, Morse, Chu, Little and Corbato 1956) of 
these functions he w^as finally able to evaluate £2(0; s) numerically. On the other 
hand, in an admirably intricate analytic tour de force Jimbo, Miwa, Mori and Sato 
(1980) expressed the Fredholm determinant of the sine kernel as 



TTS 



Es(0;s)=exp( / ^ rfx ) (1.6) 



in terms of the sigma, or Hirota, representation of the fifth Painleve equation, namely 

n n^ 



(xcr")^ + 4(xcr'-D-)(x(7'-Cr + D-'2) = 0, Cr(x) ~ - - - ^ (x ^ 0), 

1-r T'^ 



see also Mehta (2004, Chap. 21) and Tracy and Widom (2000, Sect. 4.1). With respect 
to the numerical evaluation, the latter two authors conclude in a footnote, most 
probably by comparing to Gaudin's method: "Without the Painleve representations, 
the numerical evaluation of the Fredholm determinants is quite involved." However, 
one does not need to know more than the smooth kernel sin(7r(x — y)) / [n[x — y)) 
if self to approximate £2(6; s) with the method of this paper. For instance, the Gauss- 
Legendre rule with just m = 5 quadrature points already gives, in 0.2 ms computing 
time, 15 accurate digits of the value 

£2(0,0.1) = 0.90002 7271798259 ■ ■ ■ , 

that is, by calculating the determinant of a 5 x 5-matrix easily built from the kernel. 
Even though it is satisfying to have an alternative and simpler w^ay of calculat- 
ing already known quantities, it is far more exciting to be able to calculate quantities 
that otherwise have defeated numerical evaluations so far. For instance, the joint dis- 
tribution functions of the Airy and the Airy^ processes are given as determinants 
of systems of integral operators, see Prahofer and Spohn (2002), Johansson (2003), 
Sasamoto (2005) and Borodin et al. (2007). Even though a nonlinear partial differ- 
ential equation of third order in three variables has been found by Adler and van 
Moerbeke (2005, Eq. (4.12)) for the logarithm of the joint distribution function of 
the Airy process at tw^o different times, this masterful analytic result is probably of 
next to no numerical use. And in any case, no such analytic results are yet known 
for the Airyj process. How^ever, the Nystrom-type method studied in this paper can 
easily be extended to systems of integral operators. In this way, we have succeeded 
in evaluating the tw^o-point correlation functions of both stochastic processes, see 
Section 8. 
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Outline of this paper. For the proper functional analytic setting, in Section 2 we 
review some basic facts about trace class and Hilbert-Schmidt operators. In Section 3 
we review the concept of the determinant det(I + A) for trace class operators A and 
its relation to the Fredholm determinant. In Section 4 w^e study perturbation bounds 
implying that numerical calculations of determinants can only be expected to be ac- 
curate with respect to absolute errors in general. In Section 5 w^e use the functional 
analytic formulation of the problem to obtain convergence estimates for projection 
methods of Galerkin and Ritz-Galerkin t5rpe. The convergence rate is shown to de- 
pend on a proper balance between the decay of the singular values of the operator 
and the growth of bounds on the derivatives of the corresponding singular func- 
tions. This is in sharp contrast with Section 6, where we study the convergence of 
the Nystrom-type method (1.5) by directly addressing the original definition of the 
Fredholm determinant. Here, only the smoothness properties of the kernel enter the 
convergence estimates. It turns out that, for kernels of low regularity, the order of 
convergence of the Nystrom-t5rpe method can be even higher than that of a Ritz- 
Galerkin method. In Section 7 we give examples for the exponential convergence 
rates enjoyed by analytic kernels. To this end w^e discuss the details of the numerical 
evaluation of the determinants of the sine and Airy kernels, which express the prob- 
ability distributions £2(0; s) and F2(s) (the Tracy-Widom distribution) of random 
matrix theory. Finally, after extending the Nystrom-type method to systems of inte- 
gral operators we report in Section 8 on the numerical evaluation of the two-point 
correlation functions of the Airy and Airy^ processes. 

2. Trace Class and Hilbert-Schmidt Operators. We begin by recalling some 
basic material about the spectral theory of nonselfadjoint compact operators, which 
can be found, e.g., in Gohberg, Goldberg and Kaashoek (1990), Lax (2002) and Simon 
(2005). We consider a complex, separable Hilbert space Ti. with an inner product (■, ■ ) 
that is linear in the second factor and conjugate linear in the first. The set of bounded 
linear operators will be denoted by B{H), the compact operators by Joo{'H). The 
spectrum of a compact operator A G Jcx,{'H) has no non-zero limit point; its non- 
zero points are eigenvalues of finite algebraic multiplicity. We list these eigenvalues 

as (A,i(A))jjJ^j , counting multiplicity, where N(A) is either a finite non-negative 
integer or infinity, and order them by 

|Ai(A)|;j|A2(A)|^---. 
The positive eigenvalues 

si{A) ^S2(A) ^ ■■■ > 
of the associated positive-semidefinite, selfadjoint operator 

|A| = (A*A)i/2 

are called the singular values of A. Correspondingly, there is the Schmidt or singular- 
value representation of A, that is, the norm convergent expansion 

N(\A\) 
A= Y^ Sn{A){Unr)'Vn, (2.I) 
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where the Un and v„ are certain (not necessarily complete) orthonormal sets in Ti. 
Note that s„{A) = |A„(A)| if A is selfadjornt. In general we have Weyl's inequality 

N N 

J2\K{AW ^ J2'n{A)P (N ^ N(A), 1 ^ p < «5). (2.2) 

n=l n=l 

The Schatten-von Neumann classes of compact operators are defined as 

N{\A\) 

JpiU) = {Ae JooiH) : Y. 'n{AY < <x)} (1 ^ p < <x,) 



with the corresponding norm^ 



^N(\A\\ \ 1/P 



II^IU, = ( E 'n{A)r\ . 

The operator norm on Joo{'H) perfectly fits into this setting if w^e realize that 

l|Al|=si(A)= max s„{A) = \\A\\j^. 

«=!,.. .,N{|A|) 

There are the continuous embeddings Jp{T~t) C Jq{H) ior 1 ^ p ^ q ^ 00 with 

IIAIU ^ IIAIU . 

\] \]Jq ^^ W \]Jp 

The classes Jp{Ti-) are two-sided operator ideals in B{Ti.), that is, for A E Jp{'H) 
(1 ^ p ^ 00) and B e B{n) we have AB, BA e JpiU) with 

||AB||^,, ^ \\A\\j^,\\Bl \\BA\\j^ <: \\B\\ \\A\\j^,. (2.3) 

Of special interest to us are the trace class operators J\{'H) and the Hilbert-Schmidt 
operators J2{T~(-)- The product of two Hilbert-Schmidt operators is of trace class: 

||AB||^^ ^ ||A||^J|B||j^ {A,B e J2{n)). 

The trace of a trace class operator A is defined by 

tr(A) = Y,{u„,Aun) 

n=l 

for any orthonormal basis {u„)„. A deep theorem of Lidskii's (Simon 2005, Chap. 3) 
tells us that 

N{A) 

tr(A) = Y^ Xn{A), (2.4) 

which implies by Weyl's inequality (2.2) that 

N{A) 

|tr(A)|^ Y |A„(A)|^tr(|A|) = ||A|U^. (2.5) 



^In matrix theory these norms are not commonly used — with the exception of p = 2: j|A|| j^ is then 
the Schur or Frobenius norm of the matrix A. 
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Likewise, for a Hilbert-Schmidt operator A £ J2{'H) we have 

N(A) N(A) 

tr(A2) = Y. ^n{A)\ |tr(A2)| ^ ^ |A„(A)|2 ^ \\A\\%. (2.6) 

«=1 n=l 

Integral operators with L^-kernel. In the Hilbert space H = L^{a,b) of square- 
integrable functions on a finite interval {a,b) the Hilbert-Schmidt operators are ex- 
actly given by the integral operators with L^-kernel. That is, there is a one-to-one 
correspondence (Simon 2005, Thm. 2.11) between A e J72(W) and K e L^{{a,b)-^) 
mediated through 

rb 

Au{x) = / K{x,y)u{y)dy (2.7) 



with equality of norms H^Hj^ = II^IIl^- f^e spaces J2{'H) and L'^{{a,b)'^) are thus 
isometrically isomorphic. In particular, by (2.6) and a well known basic result on 
infinite products (Knopp 1964, p. 232), we get for such operators that 

N(A) N{A) 

J^ (1 -h A„(A)) converges (absolutely) 4=> J^ A„ (A) converges (absolutely). 
«=i 11=1 

Since the product is a natural candidate for the definition of det(I -|- A) it makes 
sense requiring A to be of trace class; for then, by (2.5), the absolute convergence of 
the sum can be guaranteed. 

Integral operators with a continuous kernel. A continuous kernel K E C{[a,b]^) is 
certainly square-integrable. Therefore, the induced integral operator (2.7) defines a 
Hilbert-Schmidt operator A on the Hilbert space Ti = L^{a,b). Moreover, other than 
for L^ kernels in general, the integral 

rb 
K{x,x) dx 



over the diagonal of {a,b)^ is now well defined and constitutes, in analogy to the 
matrix case, a "natural" candidate for the trace of the integral operator. Indeed, if an 
integral operator A with continuous kernel is of trace class, one can prove (Gohberg 
et al. 2000, Thm. 8.1) 

,-b 
tr(A) = / K{x,x)dx. (2.8) 

J a 

Unfortunately, however, just the continuity of the kernel K does not guarantee the 
induced integral operator A to be of trace class. ^ Yet, there is some encouraging 
positive experience stated by Simon (2005, p. 25): 

However, the counter-examples which prevent nice theorems from 
holding are generally rather contrived so that I have found the fol- 
lowing to be true: If an integral operator with kernel K occurs in 
some 'natural' way and j \K[x,x)\dx < 00, then the operator can 
(almost always) be proven to be trace class (although sometimes 
only after some considerable effort). 



9A counter-example was discovered by Carleman (1918), see also Gohberg et al. (2000, p. 71). 
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Nevertheless, we will state some simple criteria that often w^ork w^ell: 

1. If the continuous kernel K can be represented in the form 

K{x,y) ^ Ki{x,y)K2{z,y)dz {x,y e [a,b]) 

with Kj e L^{{a,b) x (c,d)), Kj E L^{{c,d) x {a,b)), then the induced in- 
tegral operator A is trace class on L {a,b). This is, because A can then be 
written as the product of tw^o Hilbert-Schmidt operators. 

2. If K(x,y) and dyK{x,y) are continuous on [a,b]'^, then the induced integral 
operator A is trace class on L?-{a,b). This is, because we can write A by 
partial integration in the form 

•b 

'y 



Au{x) = K[x,b) I u{y)dy— / I / dzK{x,z)dz\ u{y)dy 

Ja Ja \Jy I 



as a sum of a rank one operator and an integral operator that is trace class 
by the first criterion. In particular, integral operators with smooth kernels 
are trace class (Lax 2002, p. 345). 

3. A continuous Hermitian^" kernel K(x, y) on \a, b\ that satisfies a Holder con- 
dition in the second argument with exponent a > 1/2, namely 

|K(x,i/i) -K(x,i/2)| ^ Clyi -yiT {^rV\'Vi e VM), 

induces an integral operator A that is trace class on L^{a,b); see Gohberg 
et al. (2000, Thm. IV.8.2). 

4. If the continuous kernel K induces a selfadjoint, positive-semidefinite inte- 
gral operator A, then A is trace class (Gohberg et al. 2000, Thm. IV.8.3). The 
hjrpothesis on A is fulfilled for positive-semidefinite kernels K, that is, if 

n 

E ZjZi,K{xj,Xk) ^0 (2.9) 

for any Xi,...,x„ e {a,b),z e C" and any m e N (Simon 2005, p. 24). 

3. Definition and Properties of Fredholm and Operator Determinants. In this 
section w^e give a general operator theoretical definition of infinite dimensional de- 
terminants and study their relation to the Fredholm determinant. For a trace class 
operator A £ J\{T~i) there are several equivalent constructions that all define one 
and the same entire function 

rf(z) =det(i + zA) (zeC); 

in fact, each construction has been chosen at least once, in different places of the 
literature, as the basic definition of the operator determinant: 



^°An L^-kernel K is Hermitian if K(x,y) — K{y,x). This property is equivalent to the fact that the 
induced Hilbert-Schmidt integral operator A is selfadjoint. A* = A. 
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1. Gohberg and Krein (1969, p. 157) define the determinant by the locally uni- 
formly convergent (infinite) product 

N{A) 

det{I+zA)= Yl{l+zX„{A)), (3.1) 

n=l 

which possesses zeros exactly at z„ = —1/A„{A), counting multiplicity. 

2. Gohberg et al. (1990, p. 115) define the determinant as follows. Given any 
sequence of finite rank operators A,, with A„ — > A converging in trace class 
norm, the sequence of finite dimensional determinants^^ 

det(l + zA„r,an(^^,)) (3.2) 

(which are pol5momials in z) converges locally uniform to det(J + zA), in- 
dependently of the choice of the sequence A„. The existence of at least one 
such sequence follows from the singular value representation (2.1). 

3. Dunford and Schw^artz (1963, p. 1029) define the determinant by what is 
often called Plemelj's formula ^^ 

/ ~ (-z)" \ 

det(I + zA) = exp(trlog(I + zA)) = exp - ^ ^^ ^trA" , (3.3) 

which converges for |z| < l/|Ai(A)| and can analytically be continued as 
an entire function to all z G C 

4. Grothendieck (1956, p. 347) and Simon (1977, p. 254) define the determinant 
most elegantly with a little exterior algebra (Greub 1967). With A"(^) G 
^i(A"(^)) being the n* exterior product of A, the power series 

00 

det(I + zA)=X]^"trA"(^) (34) 

converges for all z e C. Note that tr A"(^) = E,i<--<i„ A/j(A) ■ ■ ■ A/^(A) is 
just the n* symmetric function of the eigenvalues of A. 
Proofs of the equivalence can be found in (Gohberg et al. 2000, Chap. 2) and (Simon 
2005, Chap. 3). We will make use of all of them in the course of this paper. We state 
two important properties (Simon 2005, Thm. 3.5) of the operator determinant: First 
its multiplication formula, 

det(I + A + B + AB)=det(I + B)det(I + A) {A,B e Ji{n)), (3.5) 

then the characterization of invertibility: det(I + A) 7^ if and only if the inverse 
operator (I + A)~^ exists. 



'^Gohberg et al. (2000) have later extended this idea to generally define traces and determinants on 
embedded algebras of compact operators by a continuous extension from the finite dimensional case. 
Even within this general theory the trace class operators enjoy a most unique position: it is only for them 
that the values of trace and determinant are independent of the algebra chosen for their definition. On the 
contrary, if A is Hilbert-Schmidt but not trace class, by varying the embedded algebra, the values of the 
trace tr(A) can be given any complex number and the values of the determinant det(I + A) are either 
always zero or can be made to take any value in the set C \ {0} (Gohberg et al. 2000, Chap. VII). 

'^Plemelj (1904, Eq. (62)) had given a corresponding form of the Fredholm determinant for integral 
operators. However, it can already be found in Fredholm (1903, p. 384). 
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The matrix case. In Section 6 we will study the convergence of finite dimensional 
determinants to operator determinants in terms of the pow^er series (3.4). Therefore, 
we give this series a more common look and feel for the case of a matrix A G c^^™. 
By evaluating the traces with respect to a Schur basis of A one gets 

trA"(^)= E det(A,„,,)^^^^i = - ^ det(A,„,)^,^^„ 
ii<---<i„ ■ /I,. ..,;„=! 

that is, the sum of all n x n principal minors of A. The yields the von Koch (1892) 
form of the matrix determinant 

det(I + zA)=^- Y. det(A,p,,,);,,^i (A e C'«><'"). (3.6) 

n=0 "■ ii,.. .,/„=! 

(In fact, the series must terminate at n = m since det(J + zA) is a pol5n-iomial of de- 
gree ;« in z.) A more elementary proof of this classical formula, by a Taylor expansion 
of the polynomial det(I + zA), can be found, e.g., in Meyer (2000, p. 495). 

The Fredholm determinant for integral operators with continuous kernel. Suppose that 
the continuous kernel K E C{[a,b]^) induces an integral operator A that is trace class 
on the Hilbert space H = L^{a, b). Then, the traces of A" (^) evaluate to (Simon 2005, 
Thm. 3.10) 

trA"(A)=l/ det{K{t^,t^)r ,dt,---dtn (n = 0,1,2,. . .). 

ni J (a,h)" 

The power series representation (3.4) of the operator determinant is therefore exactly 
Fredholm's expression (1.2), that is, 

00 n ,■ 

det(I + zA) = ^ - / det(X(fp, f,)); 1 dti ■ ■ ■ dtn. (3.7) 

The similarity with von Koch's formula (3.6) is striking and, in fact, it was just an 
analogy in form that had led Fredholm to conceive his expression for the deter- 
minant. It is important to note, how^ever, that the right hand side of (3.7) perfectly 
makes sense for any continuous kernel, independently of whether the corresponding 
integral operator is trace class or not. 

The regularized determinant for Hilbert-Schmidt operators. For a general Hilbert- 
Schmidt operator we only know the convergence of X^„ A(A)-^ but not of X]n ^n{-^)- 
Therefore, the product (3.1), which is meant to define det(I + zA), is not known 
to converge in general. Instead, Hilbert (1904) and Carleman (1921) introduced the 
entire function^^ 

N(A) 

det2(/ + zA)= n(l + zA„(A))e-^^«('^) {A e JiiU)), 

n=\ 



^3In fact, using such exponential "convergence factors" is a classical technique in complex analysis to 
construct, by means of infinite products, entire functions from their intended sets of zeros, see Ablowitz 
and Fokas (2003, Sect. 3.6). A famous example is 



• "fl(. + -l. 



-z/ti 



r(2) ^\\ n, 

which corresponds to the eigenvalues A„(A) — 1/n of a Hilbert-Schmidt operator that is not trace class. 
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which also has the property to possess zeros exactly at Zj, = — 1/A„(A), counting 
multiplicity. Plemelj's formula (3.3) becomes (Gohberg et al. 2000, p. 167) 

det2(I + zA) = exp j - f; t^trA" J 

for \z\ < l/|Ai(A)|, which perfectly makes sense since A^, fi? , . . .are trace class for 
A being Hilbert-Schmidt.^'^ Note that for trace class operators we have 

det(I + zA) = det2(I + zA)exp(z-trA) (A (z Ji{H)). 

For integral operators A of the form (2.7) with a continuous kernel KonTC = L^{a, b) 
the Fredholm determinant (1.2) is related to the Hilbert-Carleman determinant by 
the equation (Hochstadt 1973, p. 250) 

rb 

d{z) = det2(-f + zA) exp(z / K{x, x) dx) 



in general, even if A is not of trace class. It is important to note, though, that if 

A ^ J^iCH) with such a kernel, we have J^ K{x,x) dx ^ tr(A) simply because tr(A) 
is not well defined by (2.4) anymore then. 

4. Perturbation Bounds. In studying the conditioning of operator and matrix 
determinants we rely on the fundamental perturbation estimate 

|det(I + A) -det(I + B)| ^ \\A - B\\j^exp {1 + max{\\A\\j^,\\B\\jJ) (4.1) 

for trace class operators, which can beautifully be proven by means of complex anal- 
ysis (Simon 2005, p. 45). This estimate can be put to the form 

|det(I+(A + £))-det(7 + A)| ^ e^+H^H-^i ■ ||£|| ^^ + 0(||£||2j, (4.2) 

showing that the condition number k^^s of the operator determinant det(I + A), with 
respect to absolute errors measured in trace class norm, is bounded by 



'^abs 



<;,l+l|A||jl. 



This bound can considerably be improved for certain selfadjoint operators that will 
play an important role in Section 7. 



'^Equivalently one can define (Simon 2005, p. 50) the regularized determinant by 

deiiil + zA) = det(I + {(I + zA)e-''^ - I)) (A e JiiH)). 

This is because one can then show (I + zA)e^^^ — 7 £ Ji{'H). For integral operators A on L^[a,b) with 
an L^ kernel K, Hilbert (1904, p. 82) had found the equivalent expression 

K(ti,t2) ■■■ K(h,t„] 

VI 

det2{/ + 2A)= Y,- / 



K(h,h) ■■■ K{t2,tn) 



K{t„,h) K(t,„t2) ■■■ 



dti ■ ■ ■ dtfj. 



simply replacing the problematic "diagonal" entries K{tj,tj) in Fredholm's determinant (3.7) by zero. 
Simon (2005, Thm. 9.4) gives an elegant proof of Hilbert's formula. 



ON THE NUMERICAL EVALUATION OF FREDHOLM DETERMINANTS 13 



Lemma 4.1. Let A e j7i(W) be selfadjoint, positive-semidefinite with Ai(A) < 1. If 



|E||jj < ||(I-A)-i||-i then 



|det(I-(A + £))-det(I-A)|^ \\E\\j^. (4.3) 

That is, the condition number k„j,s of the determinant det(I — A), with respect to absolute 
errors measured in trace class norm, is bounded by K^bs ^ 1- 

Proof. Because of 1 > Ai(A) ^ A2(A) ^ ■ ■ ■ ^ there exists the inverse operator 
(I — A)^^. The product formula (3.1) implies det(I — A) > 0, the multiplicativity (3.5) 
of the determinant gives 

det(I - (A + E)) = det(I - A) det(I - (I - A)-1e). 

Upon applying Plemelj's formula (3.3) and the estimates (2.3) and (2.5) we get 

{{i-A)-^Er 



|logdet(I-(I-A)-i£) 



tME 



.«=! 



^ E = log 



„-i " v^-iia-^)-iii-ii£iUi 

if 11(1 — A)~^|| ■ \\E\\j^ < 1. Hence, exponentiation yields 

1 - 11(1 - A)-^\\ ■ \\E\\j^ <; det(I - (I - A)-1e) 

-■- 111-' ^j II II'^IIji 
that is 



|det(I-(A + E))-det(I-A)| ^ det(I - A) • \\{I-A)-'^\\ ■ \\E\ 
Now, by the spectral theorem for bounded selfadjoint operators we have 

1 N(A) J ^ 



Jl- 



K'-^'-ii^r^ATW^ni 



1 - Ai(A) ^ „^i^ 1 - An{A) det(J - A) 

and therefore det(I — A) ■ 1|(I — A)^^|| ^ 1, which finally proves the assertion. D 

Thus, for the operators that satisfy the assumptions of this lemma the determi- 
nant is a really well conditioned quantity — with respect to absolute errors, like the 
eigenvalues of a Hermitian matrix (Golub and Van Loan 1996, p. 396). 

Implications on the accuracy of numerical methods. The Nystrom-type method of 
Section 6 requires the calculation of the determinant det(i + A) of some matrix A G 
C"^™. In the presence of roundoff errors, a backw^ard stable method like Gaussian 
elimination with partial pivoting (Higham 2002, Sect. 14.6) gives a result that is exact 
for some matrix A = A + E with 

|£;,)cl ^ f^l^/Vcl {i,k = l,...,m) (4.4) 
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where e is a small multiple of the unit roundoff error of the floating-point arithmetic 
used. We now use the perturbation bounds of this section to estimate the resulting 
error in the value of determinant. Since the trace class norm is not a monotone matrix 
norm (Higham 2002, Def. 6.1), we cannot make direct use of the componentwise 
estimate (4.4). Instead, we majorize the trace class norm oi my. m matrices A by the 
Hilbert-Schmidt (Frobenius) norm, which is monotone, using 

^1/2 

\A\\j, ^^\\A\\j^, \\A\\j^= ( Y. 

Thus, the general perturbation bound (4.2) yields the following a priori estimate of 
the roundoff error affecting the value of the determinant: 

I det(i + A) - det(I + A)| ^ Vm\\A\\j^ exp (l +\\A\\jJ-e + 0{e^). 

If the matrix A satisfies the assumptions of Lemma 4.1, the perturbation bound (4.3) 
gives the even sharper estimate 

I det(I -A)- det(I - A)| ^ V^\\A\\j^ ■ e. (4.5) 

Therefore, if det(I — A) -C ||A|| j^/ w^e have to be prepared that we probably cannot 
compute the determinant to the full precision given by the computer arithmetic used. 
Some digits will be lost. A conservative estimate stemming from (4.5) predicts the 
loss of at most 



^°§io [ det(I - A) ) 

decimal places. For instance, this will affect the tails of the probability distributions 
to be calculated in Section 7. 

5. Projection Methods. The general idea (3.2) of defining the infinite dimen- 
sional determinant det(7 + A) for a trace class operator A by a continuous exten- 
sion from the finite dimensional case immediately leads to the concept of a projec- 
tion method of Galerkin type. We consider a sequence of wi-dimensional subspaces 
Vm C Ti. together with their corresponding orthonormal projections 

1 til '■ it > Vjfi . 

The Galerkin projection PmAP„i of the trace class operator A is of finite rank. Given 
an orthonormal basis (pi,. . . ,(p„i of V,„, its determinant can be effectively calculated 
as the finite dimensional expression 

det(I + Z PmAPrn) = det {I + Z P,„AP,n \v,„ ) = det {Sij + z {(Pi, A(Pj)y.'.^^ (5.1) 

if the matrix elements {(pi,A(pj) are numerically accessible. 

Because of \\Pm\\ ^ 1/ and thus ||PmAP„,||jj ^ II^IIji ^ 1/ the perturbation 
bound (4.1) gives the simple error estimate 

|det(I + zP,„AP„)-det(I + zA)| ^ ||P™AP,„ - A||j^ ■ |z| e^+l'l-H^H'?!. (5.2) 

For the method to be convergent we therefore have to show that PmAPm ^ A in 
trace class norm. By a general result about the approximation of trace class operators 
(Gohberg et al. 1990, Thm. 4.3) all we need to know is that P„, converges pointwise^^ 



'5In infinite dimensions P„, cannot converge in norm since the identity operator is not compact. 
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to the identity operator I. This pointwise convergence is obviously equivalent to the 
consistency of the family of subspaces V,,,, that is, 

00 

\^ V„, is a dense subspace of H. (5.3) 

ni=l 

In summary, we have proven the following theorem. 

Theorem 5.1. Let A be a trace class operator. If the sequence of subspaces satisfies the 
consistency condition (5.3), the corresponding Galerkin approximation (5.1) of the operator 
determinant converges, 

dei{l + zPrnAPm) ^det(I + zA) {m ^00), 

uniformly for hounded z. 

A quantitative estimate of the error, that is, in view of (5.2), of the projection 
error \\PmAP„, — A\\j^ in trace class norm, can be based on the singular value rep- 
resentation (2.1) of A and its finite-rank truncation Aj^: (We assume that A is non- 
degenerate, that is, N(|A|) = 00; since otherwise we could simplify the following by 
putting An = A.) 

CO N 

A = Y^Sn{A){Un,-)Vn, A^ = ^ S„(A) (w„, ■ )37„. 

n=\ n=X 

We obtain, by using ||P„i || ^ 1 once more, 

\\PmAPm - A\\j^ <i \\PmAPm - PmA^PrnWa^ + \\PmAi^Pm - Af^\\j^ + l|Ajv - A\\j^ 

^ 2||An - A\\j^ + WPmA^Pm - AnII jj 

00 N 

^2 Yli Sn{A) + Y^Sn{A)\\{PmUn,-)PmV„- {Un,-)Vn\\j-^ 
n=N+l n=l 

^2 Y^ Sn{A) + Y^S„{A){\\Un-PmUn\\ + \\v„-Pm'V„\\). (54) 

fj=N+l »=1 

There are two competing effects that contribute to making this error bound small: 
First, there is the convergence of the truncated series of singular values, 

CO 

Y, s„(A)^0 (N^oj), 
which is independent of m. Second, there is, ior fixed N, the collective approximation 

P^Un -^ U„, PmV„ -^Vn {m ^ oo) 

of the first N singular functions u„, v„ {n = 1, . . . , N). For instance, given e > 0, w^e 
can first choose N large enough to push the first error term in (5.4) below e/2. After 
fixing such an N, the second error term can be pushed below e/2 for m large enough. 
This way we have proven Theorem 5.1 once more. However, in general the conver- 
gence of the second term might considerably slow down for growing N. Therefore, 
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a good quantitative bound requires a carefully balanced choice of N (depending on 
m), which in turn requires some detailed knowledge about the decay of the singular 
values on the one hand and of the growth of the derivatives of the singular functions 
on the other hand (see the example at the end of this section). While some general re- 
sults are available in the literature for the singular values — e.g., for integral operators 
A induced by a kernel K E C'^~^'^{[a, b]^) the bound 

s„(A) = 0(n-*^-z) (n^oo) (5.5) 

obtained by Smithies (1937, Thm. 12) — the results are sparse for the singular func- 
tions (Fenyo and Stolle 1982-1984, §8.10). Since the quadrature method in Section 6 
does not require any such knowledge, we refrain from stating a general result and 
content ourselves with the case that there is no projection error in the singular func- 
tions; that is, w^e consider projection methods of Ritz-Galerkin type for selfadjoint 
operators. 

Theorem 5.2. Let A be a selfadjoint integral operator that is induced by a continuous 
Hermitian kernel K and that is trace class on the Hilbert space H = L^{a, b). Assume that A 
is not of finite rank and let {u„) be an orthonormal basis of eigenfunctions of A. We consider 
the associated Ritz projection ?„, that is, the orthonormal projection 



Pm-.H^ V„i = spanJMi, ...,«„}• 



Note that in this case 



det{I + zP„,APm) = Yl{l+z\„{A)). 

IfK E C^~^'^{[a, b]^), then there holds the error estimate (^.2) with 
\\P„,AP„ -A\\j^= o(otz-'^) {m -^ 00). 

If K is bounded analytic on £p x £p (with the ellipse £p defined in Theorem A. 2), then the 
error estimate improves to 

\\PmAP„, - A\\j^ = 0(p-'"(i-^)/4) (m ^ 00), 

for any fixed choice ofe>0. 

Proof. With the spectral decompositions 

00 m 

A= ^ A„(A)(m,„-)m„, P„,APm = A„, = Y^\„{A){u„,-)u„, 

n=\ n=\ 

at hand the bound (5.4) simplifies, for N = m, to 

00 

\\P,nAPm-A\\j^= ^ |A„(A)|. 
n=m+l 

Now, by some results of Hille and Tamarkin (1931, Thm. 7.2 and 10.1), we have, for 
K E C'^^^'^{[a,b]^), the eigenvalue decay 

A„(A)=o(n-'^-z) («^oo) (5.6) 



ON THE NUMERICAL EVALUATION OF FREDHOLM DETERMINANTS 17 

(which is just slightly stronger than Smithies' singular value bound {3.3)) and, for K 
bounded analytic on £p x £p, 

which proves both assertions. D 

However, for kernels of low regularity, by taking into account the specifics of 
a particular example one often gets better results than stated in this theorem. (An 
example with an analytic kernel, enjoying the excellent convergence rates of the 
second part this theorem, can be found in Section 7.) 

An example: Poisson's equation. For a given / e L^(0, 1) the Poisson equation 

-u"{x)=f{x), m(0) = m(1)=0, 

with Dirichlet boundary conditions is solved (Hochstadt 1973, p. 5) by the applica- 
tion of the integral operator A, 

u{x) = Af{x) = f^ K{x,y)f{y) dy, (5.7) 

which is induced by the Green's kernel 

H^'V) = { ,, , , . (5-8) 

yy[l — x) otherwise. 

Since K is Lipschitz continuous, Hermitian, and positive definite, w^e know from the 
results of Section 2 that A is a selfadjoint trace class operator onTi. = L^(0, 1). The 
eigenvalues and normalized eigenfunctions of A are those of the Poisson equation 
which are known to be 

An(A) = „ „ , u n( x) ^ V2sm( nnx) (n = l,2,...). 

Note that the actual decay of the eigenvalues is better than the general Hille-Ta- 
markin bound (5.6) w^hich would, because of K G C'^'^([0, 1]^), just give A„(A) = 
Q^j^-3/2^ The trace formulas (2.4) and (2.8) reproduce a classical formula of Euler's, 
namely 

tr(A) = / K{x,x)dx = I x{l—x)dx 



^^ n^n^ Jo ' Jo 6' 

whereas, by (3.1) and the product representation of the sine function, the FredhoLm 
determinant explicitly evaluates to the entire function 

sin(Vz) 



^e.a--) = n(>-s^)^^. <=■') 



The sharper perturbation bound of Lemma 4.1 applies and we get, for each finite 
dimensional subspace V„j c Ti. and the corresponding orthonormal projection P^ '■ 
H -^ Vm, the error estimate 

I det(I - Pn,AP„:) - det(I - A) I ^ \\PmAP„, -A\\j^. (5.10) 

Now, w^e study two particular families of subspaces. 



i8 
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Fig. 1. Convergence of Ritz-Galerkin (circles) and Galerkin (dots) for approximating the Fredholm determinant 
of the integral operator induced by Green's kernel ofPoisson's equation. The solid line shows the upper hound l/n^m 
of Ritz-Galerkin as given in (;.ii). 



Trigonometric polynomials. Here, we consider the subspaces 

Vm = span{sin(n7r-) : n = 1,. . . ,m} = span{M„ : n = 1, . . . ,m}, 

which are exactly those spanned by the eigenf unctions of A. In this case, the projec- 
tion method is of Ritz-Galerkin type; the estimates become particularly simple since 
w^e have the explicit spectral decomposition 

oo 

PrtiAPjn — A ^ 2_^ An(^A)\U„, ■ )Ufi 

n=m+l 

of the error operator. Hence, the error bound (5.10) directly evaluates to 

\det{I-P,nAP„,)-det{I-A)\ 



^ \\Pm^Pm ~ ^\\ji 



00 1 "^ 1 1 

n=m+l n=m+l 



(5.11) 



Figure 1 shows that this upper bound overestimates the error in the Fredholm deter- 
minant by just about 20%. 

Algebraic polynomials. Here, we consider the subspaces of algebraic polynomials 
of order m, that is, 

Vm = {m G L^(0, 1) : M is a polynomial of degree at most m — 1}. 

We apply the bound given in (5.4) and obtain (keeping in mind that A is self adjoint) 

00 N 

\\PmAP„, -A\\j^^2 Y^ A„(A) + 2 ^ \„{A) \\u„ - P„,m„|| 

with a truncation index N yet to be skilfully chosen. As before in (5.11), the first term 
of this error bound can be estimated by 2/tz^N. For the second term we recall that 
the projection error || m„ — P,nU„ \\ is, in fact, just the error of pol5momial best approxi- 
mation of the eigenfunction u„ with respect to the L-^-norm. A standard Jackson-t5rpe 
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inequality (DeVore and Lorentz 1993, p. 219) from approximation theory teaches us 
that 



^n ^ m^n 



< 



^k ii„W 



m'' 



Ck 



(Tin) 



jjj/i 



where c^- denotes a constant that depends on the smoothness level k. A fixed eigen- 
function u„ (being an entire function in fact) is therefore approximated beyond every 
algebraic order in the dimension m; but with increasingly larger constants for higher 
"wave numbers" n. We thus get, with some further constant c^- depending onk ^ 2, 



\1 m-n-l m 



A\ 



^h^-'A^ 



N' 



*-i 



{k-\)mi' 



We now balance the two error terms by minimizing this bound: the optimal trunca- 
tion index N turns out to be exactly N = m, which finally yields the estimate 



det(I - PniAPm) - det(I - A)\ ^ \\PmAPm - All . ^ 



C/c 



m 



Thus, in contrast to the approximation of a single eigenfunction, for the Fredhobn 
determinant the order of the convergence estimate does finally not depend on the 
choice of k anymore; we obtain the same 0{ni^^) behavior as for the Ritz-Galerkin 
method. In fact, a concrete numerical calculation^ shows that this error estimate 
really reflects the actual order of the error decay of the Galerkrn method, see Figure 1. 

Remark. The analysis of this section has shown that the error decay of the pro- 
jection methods is essentially determined by the decay 



k=m+\ 







of the singular values of A, which in turn is related to the smoothness of the kernel 
K of the integral operator A. In the next section, the error analysis of Nystrom-type 
quadrature methods will relate in a much more direct fashion to the smoothness of 



'''By (5.1) all we need to know for implementing the Galerkin method are the matrix elements {(pi, A(pj) 
for the normalized orthogonal polynomials cp,, (that is, properly rescaled Legendre polynomials) on the 
interval [0,1]. A somewhat lengthy but straightforward calculation shows that these elements are given 
by 

\ 



i{<Pi,A<Pj))i,j = 



r- 





bo 







1 

60 





bl 


bo 





fll 






h 




H2 



■■/ 



with the coefficients 



2(2n + l)(2n + 5)' 



bn = 



4(2« + 3) V(2n + l)(2n + 5) 
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the kernel K, giving even much better error bounds, a priori and in actual numerical 
computations. For instance, the Green's kernel (5.7) of low regularity can be treated 
by an m-dimensional approximation of the determinant with an actual convergence 
rate of 0{m^^) instead of 0{m^^) as for the projection methods. Moreover, these 
methods are much simpler and straightforw^ardly implemented. 

6. Quadrature Methods. In this section we directly approximate the Fredholm 
determinant (1.2) using the Nystrom-type method (1.4) that we have motivated at 
length in Section 1. We assume throughout that the kernel K is a continuous function 
on [fl, f']^. The notation simplifies considerably by introducing the n-dimensional 
ftmctions X„ defined by 

Mh t„)=det(X(ip,f<,))j;^^^j. 

Their properties are given in Lemma A.4 of the appendix. We then write the Fred- 
holm determinant shortly as 

d{z) = \+ ^ — / K„{h,...,t„)dh--- dt„. 
„^i n\ J[a,b]" 

For a given quadrature formula 

Q(/) = EV(^/) ~ I' f{x)dx 



III I'D 

Y^Wjf{Xj) w / /(x), 

/=1 •''' 



we define the associated Nystrom-tjrpe approximation of d(z) by the expression 

dg (z) = det {3ij + z WiK{xi, Xj))'".^^ . (6.1) 

The key to error estimates and a convergence proof is the observation that we can 
rewrite dQ{z) in a form that closely resembles the Fredholm determinant. Namely, by 
using the von Koch form (3.6) of matrix determinants, the multi-linearity of minors, 
and by introducing the n-dimensional product rule (A. 3) associated with Q (see the 
appendix), we get 

^ Z" ^ / \n 

n=\ ■ j^,...,j„=l f''l 

^ Z" ^ I \" 

= 1 + > — - > Wj---Wj det I K(xj ,Xj ) I 

^,nl.^'^ '" V ^ 'P' ''''Jp,q=l 

00 n m 

= 1 + E ^ E ^;i ■ ■ ■ ^in Kni^h ^/») 

«=1 ■ h,-,in=l 
00 n 

-i+E^Q"(^")- 

n=l 

Thus, alternatively to the motivation given in the introductory Section 1, we could 
have introduced the Nystr6m-t5rpe method by approximating each of the multi- 
dimensional integrals in the power series defining the Fredholm determinant with a 
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dQ{z)-d{z)=j:-iQ"{K„ 



n=l 



Kn{ti,...,tn)dti ■ ■ ■ dt„] , (6.2) 



that is, by the exponentially generating function of the quadrature errors for the 
functions K„. The following theorem generalizes a result that Hilbert (1904, p. 59) 
had proven for a specific class of quadrature formulae, namely, the rectangular rule. 

Theorem 6.1. If a family Q„i of quadrature rules converges for continuous functions, 
then the corresponding Nystrom-type approximation of the Fredholm determinant converges, 

dQ„,{z)^d{z) {m-,00), 

uniformly for bounded z. 

Proof. Let z be bounded by M and choose any e > 0. We split the series (6.2) at 
an index N yet to be chosen, getting 

N T^n p 

\dQ„,{z)-d{z)\^J2— Ql{K„)- Kn{ti,...,tn)dt^---dtn 

+ E ^ QmiKn)- Kn{ti t,,) dt^ ■ ■ ■ dtn 



n=N+l 



Let A be the stability bound of the convergent family Q,„ of quadrature rules (see 
Theorem A.i of the appendix) and put Ai = n\ax{A,b — a). Then, by Lemma A.4, 
the second part of the splitting can be bounded by 



M" 



n=N+l 



n\ 



Ql{Kn)- I K„{ti,...,tn)dti--- dt„ 

'a,b]" 



^ E 



n=N+l 



M" / r 

-^{\Q"m{Kn)\ + \ / Kn{t^,...,t„)dt^---dt„\ 

n\ \ J[a,bY 



M" 



^ E ^(A" + (2'-«)")||X„lk»^2 Y. 



.nil 



n=N+l 



n=N+l 



n\ 



-(MAi||K||i.c»)". 



The last series converges by Lemma A. 5 and the bound can, therefore, be pushed 
below e/2 by choosing N large enough. After fixing such an N, w^e can certainly 
also push the first part of the splitting, that is. 



^ M" 



n=l 



n\ 



Ql{K„) 



a,b]" 



K„{t-[,. . . ,t„) dt-[- ■ ■ dt„ 



below e/2, now for m large enough, say m ^ wzq, using the convergence of the 
product rules QJJ, induced by Q,,, (see Theorem A. 3). In summary we get 

\dQ,„{z) - d{z)\ ^ e 

for all |z| ^ M and m ^ ntQ, which proves the asserted convergence of the Nystrom- 
type method. D 
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If the kernel K enjoys, additionally, some smoothness, we can prove error esti- 
mates that exhibit, essentially, the same rates of convergence as for the quadrature 
of the sections x i— > K{x,y) and y h^ K{x,y). 

Theorem 6.2. If K & C^~^'^{[a,b]^), then for each quadrature rule Q of order v ^ k 
with positive weights there holds the error estimate 

\dQ{z) - d{z)\ ^ c,2\b - a) V^OdzKf, - ,)||K||,) , 

where cj- fs the constant (depending only on k) from Theorem A.2, and \\K\\)^ and O are the 
norm and function defined in (A. 6) and (A. 10), respectively. 

If K is bounded analytic on £p x £p (with the ellipse £p defined in Theorem A.z), then 
for each quadrature rule Q of order v with positive weights there holds the error estimate 

\dQ{z)-d{z)\^^^^<b{\z\{b-a)\\K\l^^s,.e,)). 

Proof. By Theorem A. 3 and Lemma A.4 we can estimate the error (6.2) in both 
cases in the form 

00 {«+2)/2 

\dQ{z) - d{z)\ ^ cc E —^^ {\zW = «*(|z|/5) ; 

n=\ 

with the particular values a = cj. 2^{b — a)^v^^ and [i = {b — a) \\K\\j^ in the first case, 
and a = ip^" /{I — p^^) and /5 = {b — a) l|i^l|L~(£px£p) in the second case. This proves 
both assertions. D 

An example with an analytic kernel, enjoying the excellent convergence rates of 
the second part this theorem, can be found in Section 7. 

Note that Theorem 6.2 is based on a general result (Theorem A.2) about quadra- 
ture errors that stems from the convergence rates of pol5momial best approximation. 
There are cases (typically of low regularity), however, for which certain quadrature 
formulae enjoy convergence rates that are actually better than best approximation. 
The Nystrom-type method inherits this behavior; one w^ould just have to repeat the 
proof of Theorem 6.2 then. We refrain from stating a general theorem, since this 
w^ould involve bounds on the highest derivatives rnvolvrng weights^^ that take into 
account the boundary of the interval [a, b] . Instead, we content ourselves with the 
detailed discussion of a particular example. 

An example: Poisson's equation. We revisit the example of Section 5, that is the in- 
tegral operator (5.7) belonging to the Green's kernel K (defined in (5.8)) of Poisson's 
equation. Recall from (5.9) that 

rf(-l) = det(I-A) =sin(l). 

The kernel K is just Lipschitz continuous, that is, K e C'''^([0, 1]^). If we apply the 
Nystrom-t5^e method with the wi-point Gauss-Legendre (order v = 2m) or the 
Curtis-Clenshaw (order v = m) formulae as the underlying quadrature rule Qm, 
Theorem 6.2 proves an error bound of the form 

dQ„,{-l)-d{-l) = 0{m-^), 



'7For the interval [—1, 1] this weight would be (1 — x^Y''^, see Davis and Rabinowitz (1984, §4-8. 1) 
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Fig. 2. Convergence of the Nystrom-type method for approximating the Fredholm determinant of the integral 
operator induced by Green's kernel (^.8) ofPoisson's equation; the underlying quadrature rules Q„ are the m-point 
Gauss-Legendre (dots) and Clenshaw-Curtis (circles) rules. Note that, in accordance with Trefethen (2008), both 
behave essentially the same. The solid line shows the function l/25m^, just to indicate the rate of convergence. For 
comparison we have included the results of the Ritz-Galerkin method (stars) from Figure 1. 



which superficially indicates the same convergence rate as for the m-dimensional 
Galerkin methods of Section 5. However, the actual numerical computation shown 
in Figure 2 exhibits the far better convergence rate of 0{m ). This deviation can be 
understood in detail as follows: 

On the one hand, by inverse theorems of approximation theory (DeVore and 
Lorentz 1993, p. 220), valid for proper subintervals of [a,h\, the polynomial best ap- 
proximation (of degree m) of sections of the Green's kernel K cannot give a better 
rate than 0{m~^); since otherwise those sections could not show jumps in the first 
derivative. Given the line of arguments leading from pol5n-iomial best approxima- 
tion to Theorem 6.2, the error estimate of 0{m^^) w^as therefore the best that could 
be established this way. 

On the other hand, the sections of the Green's kernel look like piecewise linear 
hat functions. Therefore, the coefficients am oi their Chebyshev expansions decay 
as 0{m^-^) (Davis and Rabinowitz 1984, Eq. (4.8.1.3)). Given this decay rate, one can 
then prove — see, for Gauss-Legendre, Davis and Rabinowitz (1984, Eq. (4.8.1.7)) and, 
for Clenshaw-Curtis, Riess and Johnson (1972, Eq. (9)) — that the quadrature error is 
of rate 0(ot^^), too. Now, one can lift this estimate to the Nystrom-like method 
essentially as in Theorem 6.2; thus proving in fact that 

dQ,„{-l)-d{-l) = 0{m-^), 

as numerically observed. 

Remark. This "superconvergence" property of certain quadrature rules, as op- 
posed to best approximation, for kernels with jumps in a higher derivative is there- 
fore also the deeper reason that the Nystrom-type method then outperforms the 
projection methods of Section 5 (see Figure 2): Best approximation, by direct (Jack- 
son) and inverse (Bernstein) theorems of approximation theory, is strongly tied with 
the regularity of K. And this, in turn, is tied to the decay of the singular values of 
the induced integral operator A, which governs the convergence rates of projection 
methods. 
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Fig. 3. The probability £2(0; s) that an interval of length s does not contain, in the bulk scaling limit of 
level spacing 1, an eigenvalue of the Gaussian unitary ensemble (GUE). The result shown was calculated with the 
NystrSm-like method based on Gauss-Legendre with m = 30; and cross-checked against the asymptotic expansion 
logE2(0;s) = -7r2s2/8-log(s)/4 + log(2)/3-log(7r)/4 + 3^'(-l) +O(s-i)/ors^ cx> (Deift et al. 1997). 



A note on implementation. If the quadrature weights are positive (which in view 
of Theorem A.i is anjrway the better choice), as is the case for Gauss-Legendre and 
Clenshaw-Curtis, w^e recommend to implement the Nystrom-type method (6.1) in 
the equivalent, symmetric form 



dQ{z) = det(I + zA 



Q) 



^Q = 



w 



1/2 



-Z\. I -mL 1 f J\. 1 f 



W: 



1/2 



W = l 



(6.3) 



(Accordingly short Matlab and Mathematica code is given in the introductory Sec- 
tion 1.) The reason is that the m x wi-matrix Aq inherits some important structural 
properties from the integral operator A: 

• If A is selfadjoint, then Aq is Hermitian (see Footnote 10). 

• If A is positive semidefinite, then, by (2.9), Aq is positive semidefinite, too. 

This way, for instance, the computational cost for calculating the finite-dimensional 
determinant is cut to half, if by structural inheritance / + zAq is Hermitian positive 
definite; the Cholesky decomposition can then be employed instead of Gaussian 
elimination with partial pivoting. 

7. Application to Some Entire Kernels of Random Matrix Theory. In this sec- 
tion w^e study tw^o important examples, stemming from random matrix theory, with 
entire kernels. By Theorem 6.2, the Nystrom-type method based on Gauss-Legendre 
or Curtis-Clenshaw quadrature has to exhibit exponential convergence. 

7.1. The sine kernel. The probability £2(8/ s) (shown in Figure 3) that an in- 
terval of length s does not contain, in the bulk scaling limit of level spacing 1, an 
eigenvalue of the Gaussian unitary ensemble (GUE) is given (Mehta 2004, Sect. 6.3) 
by the Fredholm determinant 

£2(0; s) =det(I-A,) 

of the integral operator As on L^(0, s) that is induced by the sine kernel K: 

sin(7T(x — y)) 
n{x-y) 



Aguix) = I K(x,y)u(y)dy, K(x,y) 
Jo 
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Note that K{x,i/) is Hermitian and entire on C x C; thus Ag is a selfadjoint operator 
of trace class on L^{0,s). (This is already much more than we would need to know 
for successfully applying and understanding the Nystrom-type method. However, 
to facilitate a comparison with the Ritz-Galerkin method, we analyze the operator 
As in some more detail.) The factorization 

K{x,y) = ^ r ^'^'"'^^ ^? = 7^ /" '^"'^'^''"^d^ (7-1) 

of the kernel implies that A, is positive definite with maximal eigenvalue Ai(As) < 1; 
since, for 7^ m e L^(0, s), we obtain 



< {u,Asu) = [ 

J — n 



''''hi{x)dx 



2 



d^ 



V2ti Jo 

-TT J —CO 



Here, in stating that the inequalities are strict, we have used the fact that the Fourier 
transform u of the function u G L^(0, s), which has compact support, is an entire 
function. Therefore, the perturbation bound of Lemma 4.1 applies and w^e obtain, 
for Ritz-Galerkin as for any Galerkin method, like in the example of Section 5, the 
basic error estimate 

I det(I - P,nAsPm) - det(I - A,)| ^ ||P,„A,P„, - A,|| j^. 

Now, Theorems 5.2 and 6.2 predict a rapid, exponentially fast convergence of the 
Ritz-Galerkin and the Nystrom-type methods: In fact, an wi-dimensional approxima- 
tion will give an error that decays like 0{e^"^), for any fixed c > since, for entire 
kernels, the parameter p > 1 can be chosen arbitrarily in these theorems. 

Details of the implementation of the Ritz-Galerkin method. There is certainly no gen- 
eral recipe on how to actually construct the Ritz-Galerkin method for a specific 
example, since one would have to know, more or less exactly, the eigenvalues of A. 
In the case of the sine kernel, however, Gaudin (1961) had succeeded in doing so. 
(See also Katz and Sarnak (1999, p. 411) and Mehta (2004, pp. 124-126).) He had 
observed that the integral operator At on L^( — 1, 1), defined by 



Afu{x)= I e'''^''yu{y)dy 



(which is, by (7.1), basically a rescaled "square-root" of Ait), is commuting with the 
selfadjoint, second-order differential operator 

Lu{x) = — ({x^ - l)u'{x)\ + K^t^x^u{x) 
with boundary conditions 

(1 - x2)m(x)|.,=±i = (1 - x2)w'(x)U=±i = 0. 

Thus, both operators share the same set of eigenf unctions Un, namely the radial 
prolate spheroidal wave functions (using the notation of Mathematica 6.0) 

Un{x) = S^^l{nt,x) (m = 0,1,2...). 
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Fig. 4. Convergence of various m-dimensional approximations of the Fredholm determinants £2(0;!) (left) 
and £2(0; 2) (right): the Nystrom-type quadrature methods based on Gauss-Legendre (dots) and Curtis-Clenshaw 
(circles), as well as Gaudin's (ig6i) Ritz-Galerkin method based on spheroidal wave functions (stars). The dashed 
line shows the amount, according to (4.5), of roundoff error due to the numerical evaluation of the finite-dimensional 
determinants; all calculations were done in IEEE double arithmetic with machine precision e = 2^^^. 



These special functions are even for n even, and odd for n odd. By plugging them 
into the integral operator At Gaudin had obtained, after evaluating at x = 0, the 
eigenvalues 



^ik{At 



M2/c(0) J -I 



U2k{y)dy, ^ik+i{At) 



int 



^2k+l 



(0) 



U2k+i{y)ydy. 



Finally, we have (starting with the index n = here) 



A„(A,) = -|A„(A,/2)|2 (n = 0,1,2,...). 
Hence, the m-dimensional Ritz-Galerkin approximation of det(I — As) is given by 

m—l 

det(I - P„,AsP«) = n (1 - M^s)). 

n=0 

While Gaudin himself had to rely on tables of the spheroidal wave functions (Stratton 
et al. 1956), w^e can use the fairly recent implementation of these special functions 
by Falloon, Abbott and Wang (2003), w^hich now comes with Mathematica 6.0. This 
w^ay, we get the following implementation: 

u[t_, n_, x_] : = SpheroidalSl[n, 0, 7rt, x] 

NIntegrate[u[t, n, y] , {y, -1, 1}] 



Ai[t_, n_?EvenQ] : = 
^[t_, n_?OddQ] := 



A[s_, n_] : = A[s, n] = — Abs Ui — , n 



u[t, n, 0] 
i TT t NIntegrate[y u[t, n, y] , {y, -1, 1}] 

u*"'"'!' [t, n, 0] 
2 



DetRlt zGalerkln 



[s_, m_] :=P](1-A[s, n] ) 



Given all this, one can understand that Jimbo et al.'s (1980) beautiful discovery of 
expressing £2(0; s) by formula (1.6) in terms of the fifth Painleve transcendent was 
generally considered to be a major break-through even for its numerical evaluation 
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Fig. 5. The Tracy-Widom distribution F2(s); that is, in the edge scaling limit, the probability of the maximal 
eigenvalue of the Gaussian unitary ensemble (GUE) being not larger than s. The result shown was calculated 
with the Nystrom-like method based on Gauss-Legendre with m = 50; and cross-checked against the asymptotic 
expansion logF2(-s) = -5^/12 - log(s)/8 + log(2)/24 + C"'(-l) + 0(5^3/2^^^,^ s ^ 00 (Deift et al. 2008). 



(Tracy and Widom 2000, Footnote 10). However, note how much less knowledge 
suffices for the application of the far more general Nystrom-type method: continuity 
of K makes it applicable, and K being entire guarantees rapid, exponentially fast 
convergence. That is all. 

An actual numerical ex-periment. Figure 4 shows the convergence (in IEEE machine 
arithmetic) of an actual calculation of the numerical values £2(0;!) and £2(0; 2). 
We observe that the Nystrom-type method based on Gauss-Legendre has an expo- 
nentially fast convergence rate comparable to the Ritz-Galerkin method. Clenshaw- 
Curtis needs a dimension m that is about twice as large as for Gauss-Legendre to 
achieve the same accuracy. This matches the fact that Clenshaw-Curtis has the or- 
der V = m, which is half the order v — 2m of Gauss-Legendre, and shows that 
the bounds of Theorem 6.2 are rather sharp with respect to v (there is no "kink" 
phenomenon here, cf. Trefethen (2008, p. 84)). The dashed line shows the amount, 
as estimated in (4.5), of roundoff error that stems from the numerical evaluation of 
the finite dimensional m x wi-determinant itself. Note that this bound is essentially 
the same for all the three methods and can easily be calculated in course of the 
numerical evaluation. We observe that this bound is explaining the actual onset of 
numerical "noise" in all the three methods reasonably w^ell. 

Remark. Note that the Nystrom-type method outperforms the Ritz-Galerkin me- 
thod by far. First, the Nystrom-type method is general, simple, and straightforw^ardly 
implemented (see the code given in Section 1); in contrast, the Ritz-Galerkin depends 
on some detailed knowledge about the eigenvalues and requires numerical access 
to the spheroidal wave functions. Second, there is no substantial gain, as compared 
to the Gauss-Legendre based method, in the convergence rate from knowing the 
eigenvalues exactly. Third, and most important, the computing time for the Ritz- 
Galerkin runs well into several minutes, w^hereas both versions of the Nystrom-type 
method require just a few milliseconds. 

7.2. The Airy kernel. The Tracy-Widom distribution F2(s) (shown in Figure 5), 
that is, in the edge scaling limit, the probability of the maximal eigenvalue of the 
Gaussian unitary ensemble (GUE) being not larger than s, is given (Mehta 2004, 
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§24.2) by the Fredholm determinant 

F2(s) = det(i - As) (7.2) 

of the integral operator As on L^(s,oo) that is induced by the Airy kernel K: 

. / ^ /■~Tx/ X ^ xj Tx/ N A.i(x)Ai'(y) -Ai(y)Ai'(x) ^ , 

Asu{x)= K{x,y)u{y)dy, K{x,y) = ^-^ ^^' „ ^^^ ^. (7-3) 

J s X — y 

Note that K is Hermitian and entire on C x C; thus As is selfadjoint. (Again, this 
is already already about all we would need to know for successfully applying and 
understanding the Nystrom-type method. How^ever, we would like to show that, as 
for the sine kernel, the strong perturbation bound of Lemma 4.1 applies to the Airy 
kernel, too.) There is the factorization (Tracy and Widom 1994, Eq. (4.5)) 

K{x,y)= / Ai(x + ^)Ai(y + ^)rf^, 
Jo 

which relates the Airy kernel with the Airy transform (Vallee and Soares 2004, §4.2) 
in a similar way as the sine kernel is related by (7.1) with the Fourier transform. This 
proves, because of the super-exponential decay of Ai(x) — > as x — > 0, that As is the 
product of two Hilbert-Schmidt operators on L^(s, 00) and therefore of trace class. 
Moreover, As is positive semi-definite with maximal eigenvalue Ai(A) ^ 1; since 
by the Parseval-Plancherel equality (Vallee and Soares 2004, Eq. (4.27)) of the Airy 
transform we obtain, for u G L^(s, 00), 



^ (m,Asm) = 



Ai(x + (^)u{x) dx 




< 



2 



d^ 



Ai(x + ^)u{x)dx 



d?=\\u\\l2. 



More refined analytic arguments, or a carefully controlled numerical approximation, 
show the strict inequality Ai(A) < 1; the perturbation bound of Lemma 4.1 applies. 
Modification of the Nystrom-type method for infinite intervals. The quadrature meth- 
ods of Section 6 are not directly applicable here, since the integral operator As is 
defined by an integral over the infinite interval (s,oo). We have the following three 
options: 

1. Using a Gauss-type quadrature formula on (s, 00) that is tailor-made for the 
super-exponential decay of the Airy function. Such formulae have recently 
been constructed by Gautschi (2002). 

2. Truncating the integral in (7.3) at some point T > s. That is, before using 
the Nystrom-type method with a quadrature formula on the finite interval 
[s, T] (for which the second part of Theorem 6.2 is then applicable, showing 
exponential convergence), we approximate the Fredholm determinant (7.2) 
by 



det(I - PtAsPt) = det i^I - As l^if^sj) 

where the orthonormal projection Pj : L^(s,oo) — > L^(s, T), Pu = u ■ X\s,T]' 
denotes the multiplication operator by the characteristic function of [s, T] 
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truncation point T 



Fig. 6. Values of the expression \\PtAsPt — As\\j^, which bounds, by (7.4), the error in F2(s) committed i 
truncating the integral in (j.^) at a point T > s. 



This way we commit an additional truncation error, which has, by passing 
through the perturbation bound of Lemma 4.1, the computable bound 



I det(I - PtAsPt) - det(I - A,)| ^ \\PtAsPt - A,\\j^ ^ 



\PtAsPt-A, 



s\\J2 



\K{x,y)\^ dxdy 



1/2 



(74) 



Figure 6 shows this bound as a function of the truncation point T. We ob- 
serve that, for the purpose of calculating (within IEEE machine arithmetic) 
F2(s) for s G [—8,2] — as shown in Figure 5 — , a truncation point at T = 16 
w^ould be more than sufficiently safe. 

3. Transforming the infinite intervals to finite ones. By using a monotone and 
smooth transformation ^s : (0,1) -^ (s, 00), defining the transformed inte- 
gral operator Ag on L^(0, 1) by 



gives the identity 

F,(s) = det(l-As 



L2(s,oo) 



det(l-As\ 



L2(0,l) 



For the super-exponentially decaying Airy kernel K we suggest the transfor- 
mation 



<^s(O = s + 10tan(7r^/2) (^e(0,l)). 



(7-5) 



Note that though K[^,rj) is a smooth function on [0,1] it possesses, as a 
function on C x C, essential singularities on the lines I, = \ or rj — 1. Hence, 
we can only apply the first part of Theorem 6.2 here, which then shows, for 
Gauss-Legendre and Clenshaw-Curtis, a super-algebraic convergence rate, 
that is, 0{m ) for arbitrarily high algebraic order k. The actual numerical 
experiments reported in Figure 7 show, in fact, even exponential conver- 
gence. 
From the general-purpose point of view, we recommend the third option. It is 
straightforw^ard and does not require any specific knowledge, or construction, as 
w^ould the first and second option. 



30 



R BORNEMANN 




•.--•*V..;.n-"sS 



40 45 50 



Fig. 7. Convergence of the m-dimensional Nystrom-type approximation — using the transformation (y.;) — of 
the Fredholm determinants F2{~2) (left) and Fzi—i) (right), based on Gauss-Legendre (dots) and Curtis-Clenshaw 
(circles). The dashed line shows the amount, according to (4.5), of roundoff error due to the numerical evaluation of 
the finite-dimensional determinants; all calculations were done in IEEE double arithmetic (e = 2^^^). 



Remarks on other numerical methods to evaluate F2(s). As for the sine kernel, there is 
a selfadjoint second-order ordinary differential operator commuting with Ag (Tracy 
and Widom 1994, p. 166). Though this has been used to derive some asymptotic 
formulas, nothing is known in terms of special functions that would enable us to 
base a Ritz-Galerkin method on it. As Mehta (2004, p. 453) puts it: "In the case of 
the Airy kernel . . . the differential equation did not receive much attention and its 
solutions are not known." 

Prior to our work of calculating F2(s) directly from its determinantal expression, 
all the published numerical calculations started with Tracy and Widom's (1994) re- 
markable discovery of expressing F2(s) in terms of the second Painleve transcendent; 
namely 



F2(s) = exp ( — / (z — s)(j(z)^ dz 



with q{z) being the Hastings-McLeod (1980) solution of Painleve II, 

q"{z) =2q{z)'^ + zq{z), q{z) ^ Ai{z) as z 



00. 



(7-6) 



Initial value methods for the numerical integration of (7.6) suffer from severe sta- 
bility problems (Prahofer and Spohn 2004). Instead, the numerically stable w^ay of 
solving (7.6) goes by considering q{z) as a connecting orbit, the other asymptotic 
state being 



q{z) 



as z 



and using numerical tw^o-point boundary value solvers (Dieng 2005). 

8. Extension to Systems of Integral Operators. We now consider an N x N 
system of integrals operators that is induced by continuous kernels K,-, G C(I/ x I A 
{i,j — 1,...,N), where the 1/ C IR denote some finite intervals. The corresponding 
system of integral equations 



Ui{x) 



N „ 

+ ^E / ^{^'y)^j{y)dy = Mx) {x e k, i,j = i,...,N) (8.1) 
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defines, with u = (mj, . . . , mjv) and / = (/i, . . ■ ,fN)f an operator equation 

u + zAu = f 
on the Hilbert space H = L'^{h) © ■ ■ ■ © L^{In)- 

8.1. The Fredholm determinant for systems. Assuming A to be trace class, let 
us express det(I + zA) in terms of the system (K/,) of kernels. To this end we show 
that the system (8.1) is equivalent to a single integral equation; an idea that, essen- 
tially, can already be found in the early w^ork of Fredholm (1903, p. 388). To simplify 
notation, we assume that the Ij- are disjoint (a simple transformation of the system of 
integral equations by a set of translations will arrange for this). We then have^^ 

N 
/c=l 

by means of the natural isometric isomorphism 

N 
(mi,...,mn) ^U= J^XkUk 
k=l 

where Xk denotes the characteristic function of the interval 4- Given this picture, the 
operator A can be viewed being the integral operator on L^(I) that is induced by the 
kernel 

N 

^x,y) = Y^ Xi{x)Kij{x,y)xj{y). 
By (3.7) we finally get (cf. Gohberg et al. (2000, Thm 6.1)) 
det{I + zA) = ^ - dei(K{t,,t,)f dh---dt„ 

CO n I- / N \ 

= E^L E Xr,{h)---Xr„{tn))det{K{t,,t,))" dH---dt„ 

n=0 "■ •'' \ii,...,!„=l / 



00 n N 



E^. E y, ... , det{K{t„t,))l^^^dH---dt„ 



n=0"' h,...,i„=l-"ii>'->'h, 
<x> n N 



S^„.E.t -./"('■'""'■''')"-"' 



dt„. 



'^The general case could be dealt with by the topological sum, or coproduct, of the intervals Ij:, 

N N 

One would then use (Johansson 2003) the natural isometric isomorphism 

N / N 

k=l \k=l 
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By eventually transforming back to the originally given, non-disjoint intervals 4, the 
last expression is the general formula that we have sought for: det(I + zA) — d{z) 
with 



oo n N 



d(-) = Lh E 






dti ■ ■ ■ dtn- 



(8.2) 



This is a perfectly well defined entire function for any system X,, of continuous 
kernels, independently of whether A is a trace class operator or not. We call it the 
Fredholm determinant of the system. 

The determinant of block matrices. In preparation of our discussion of Nystrom- 
t5^e methods for approximating (8.2) we shortly discuss the determinant of N x N- 
block matrices 



/ 



A 



A 



V 



11 



Nl 



Ain\ 



^nn/ 



eC 



MxM 



A,-,- e C'"'^™', M = mi + ---+mN. 



Starting with von Koch's formula (3.6), an argument^^ that is similar to the one that 
has led us to (8.2) yields 



~ 7« N ""h '"in , .„ 

det(I + zA)=E^ E E---Edet(A,,,V, , 

«=0 "■ /i,...,/„=l/ci=l k„=l ^ ^M-i 



(8.3) 



8.2. Quadrature methods for systems. Given a quadrature formula for each of 
the intervals I,, namely 



Qi{f) = L^Oijf{xij) 



/(x) dx, 



(8.4) 



we aim at generalizing the Nystrom-type method of Section 6. We restrict ourselves 
to the case of positive weights, if,, > 0, and generalize the method from the single 
operator case as given in (6.3) to the system case in the following form: 



A 



111 



dQ{z) = det{I + zAQ), Aq 



V 



Nl 



Mn\ 
^nnJ 



with the sub-matrices A/, defined by the entries 



{'^ij)p,q — ^"ip ^ij{^ip'^jq)'^jq 



(p = l,...,OTj, q = l,...,mj). 



(8.5) 



'9 Alternatively, we can use (3.4) and, recursively, the "binomial" formula (Greub 1967, p. 121) 



A'(^o e Vi) = (A' Vo) ® (a'"' ^1) 



of exterior algebra, which is valid for general vector spaces Vg and Vj . 
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This can be as straightforwardly implemented as in the case of a single operator. 
Now, a convergence theory can be built on a representation of the error dQ{z) — d{z) 
that is analogous to (6.2). To this end w^e simplify the notation by introducing the 
following functions on I;^ x ■ ■ ■ x I,^, 



Kk .„{h,...,tn) = det i^K,^,i,^{tp,t^)^ ^^ 

and by defining, for functions / on 1;^ x ■ ■ ■ x 7,^^ , the product quadrature formula 

/ n \ '"ii 

\k=l J 71=1 j„=l 



Thus, w^e can rewrite the Fredholm determinant (8.2) in the form 



/ f{ti,...,t„)dti ■ ■ ■ dt„. 



00 « N 



„=1 "• i^,...,i„=i-'\y-yh. 



rf(z)=l+^- Y. I, ^, ,,, Kh,..J„{h tn)dh---dtn. 



Likewise, by observing the generalized von Koch formula (8.3), we put the defini- 
tion (8.5) of dQ{z) to the form 



CO n N / n \ 

rfQ(2) = i+E^ E nQj(^. J- 

Thus, once again, the Nystrom-type method amounts for approximating each mul- 
tidimensional integral of the pow^er series of the Fredholm determinant by using 
a product quadrature rule. Given this representation. Theorem 6.2 can straightfor- 
wardly be generalized to the system case: 

Theorem 8.1. IfKjj e C'^~^'^(I; x I,), then for each set (8.4) of quadrature formulae of 
a common order v ^ k with positive weights there holds the error estimate 



JQ 



(z) - d{z) = 0{v-'') (t/^00). 



uniformly for bounded z. 

If the Kij are bounded analytic on £p{Ii) x £p{Ij) (with the ellipse £p{Ii) defined, with 
respect to f, as in Theorem A. 2), then for each set (8.4) of quadrature formulae of a common 
order v with positive weights there holds the error estimate 

dQ{z) - d{z) = 0{p-') (v^~), 

uniformly for bounded z. 

8.3. Examples from random matrix theory. Here, we apply the Nystrom-type 
method (8.5) to two 2 x 2-systems of integral operators that have recently been stud- 
ied in random matrix theory. 
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Fig. 8. Values of the two-point correlation function cov(^2(0/^2(0)) of the Airy process Aiit) (solid line). 
The dashed line shows the first term of the asymptotic expansion cov(^2 (0/^2(0)) ~ t^^ as t ^ oo. 



Two-point correlation of the Airy process. The Airy process -42(0 describes, in a 
properly rescaled limit of infinite dimension, the maximum eigenvalue of Hermi- 
tian matrix ensemble whose entries develop according to the Ornstein-Uhlenbeck 
process. This stationary stochastic process was introduced by Prahofer and Spohn 
(2002) and further studied by Johansson (2003). These authors have shown that the 
joint probability function is given by a Fredholm determinant; namely 



F{A2{t) ^ si,^2(0) ^ S2) = det I 



Ao At 



lL2(si,oo)©L2(s2,t») 



with integral operators Af that are induced by the kernel functions 

{r-CO 
j^ e-S'Ai(x + f)Ai(y + ^)df, i > 0, 

- I e^^^A\{x + ^)Ai(i/ + f ) d^, otherwise. 

7 — CO 

Of particular interest is the two-point correlation function 

COV(^2(0."42(0)) = E(^2(0^2(0)) - E(^2(0)E(^2(0)) 

d'^¥{A2{t) ^ Si, ^2(0) ^82) 



(8.6) 



(8.7) 



(8.8) 



R2 



S1S2- 



3si3s 



iot.2 



■ dsids2 — c^, 



where c^ denotes the expectation value of the Tracy-Widom distribution (7.2). We 
have calculated this correlation function for ^ f ^ 100 in steps of 0.1 to an absolute 
error of ±10^^*^, see Figure 8.^° Here are some details about the numerical procedure: 

^°A table can be obtained from the author upon request. Sasamoto (2005, Fig. 2) shows a plot (which 
differs by a scaling factor of two in both the function value and the time t) of the closely related function 



giit) = y^varM2(0--42(0))/2 = y/var(A(0)) - cov(A2{t),A2m 
-without, however, commenting on either the numerical procedure used or on the accuracy obtained. 
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• Infinite intervals of integration, such as in the definition (8.7) of the kernels 
or for the domain of the integral operators (8.6) themselves, are handled by 
a transformation to the finite interval [0,1] as in Section 7.2. 

• The kernels (8.7) are evaluated, after transformation, by a Gauss-Legendre 
quadrature. 

• The joint probability distribution (8.6) is then evaluated, after transforma- 
tion, by the Nystrom-type method of this section, based on Gauss-Legendre 
quadrature. 

• To avoid numerical differentiation, the expectation values defining the two- 
point correlation (8.8) are evaluated by truncation of the integrals, partial 
integration, and using a Gauss-Legendre quadrature once more. 

Because of analyticity, the convergence is always exponential. With parameters care- 
fully (i.e., adaptively) adjusted to deliver an absolute error of ±10^^*^, the evaluation 
of the tw^o-point correlation takes, for a single time t and using a 2 GHz PC, about 20 
minutes on average. The results were cross-checked, for small t, with the asymptotic 
expansion (Prahofer and Spohn 2002, Hagg 2007) 

cov{A2{t),A2{0)) = var(^2(0)) - |var(^2(0 - ^2(0)) 

= var(^2(0)) - t + 0{t^) {t -^ 0), 

and, for large t, with the asymptotic expansion^^ (Widom 2004, Adler and van 
Moerbeke 2005) 

cov(^2(0.^2(0)) = r^ + ct-^ + o{t-^) {t ^ «.), 

where the constant c = —3.542 ■ ■ ■ can explicitly be expressed in terms of the 
Hastings-McLeod solution (7.6) of Painleve II. 

Two-point correlation of the Airy-^ process. Sasamoto (2005) and Borodin et al. (2007) 
have introduced the Airy^ process Ai{t) for which, once again, the joint probability 
distribution can be given in terms of a Fredholm determinant; namely 

P(^l(0^Si,^l(0)^S2)=detf7- f^°^ ^Jb(si,co)©L2(«2,co)j 

with integral operators At that are now induced by the kernel functions 

fAi(x + y+^^).^(-+y)+^^V3_ exp(-(x--y)V(40) ^ ^ ^ ^^ 
Kt{x,y) = I V^t 

[ Ai(x + y + f2)gt(x+y)+2t3/3^ otherwise. 



^' Adler and van Moerbeke (2005) have derived this asymptotic expansion from the masterfully ob- 
tained result that G(f,x,y) = logP{A2{t) ^ x,A2{0) ^ y) satisfies the following nonlinear 3rd order 
PDF with certain (asymptotic) boundary conditions: 

dt[dx2 dy^ J dx^dif [ dy^ ^ dxdy dx^ ^ ^ ^ 

d^G (d'^G d^G d^G 2^ fd^G d d^Gd\fd d ,^ 



dy^dx \ dx^ dxdy dy^ ■ s J y g^.3 ^y 3^3 dx J \dx dy j 
The reader should contemplate a numerical calculation of the two-point correlation based on this PDF, 
rather than directly treating the Fredholm determinant as suggested by us. 
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0.2 0.4 0.6 0.8 1 1.2 



Fig. 9. Values of the two-point correlation function cav{Ai{t),Ai(0)) of the Airy-^ process Ai{t). 

By basically employing the same numerical procedure as for the Airy process, we 
have succeeded in calculating the two-point correlation function cov {Ai{t) , Ai{0)) 
for ^ f ^ 2.5 in steps of 0.025 to an absolute error of ±10^^*^, see Figure 9.^^ For 
a single time t the evaluation takes about 5 minutes on average (using a 2 GHz PC). 
This numerical result has been used by Bornemann, Ferrari and Prahofer (2008) as 
a strong evidence that the Airy^ process is, unlike previously conjectured, not the 
limit of the largest eigenvalue in GOE matrix diffusion. 

A. Appendices. 

A.i. Quadrature Rules. For the ease of reference, we collect in this appendix 
some classical facts about quadrature rules in one and more dimensions. 

Quadrature rules in one dimension. We consider quadrature rules of the form 



Q(/) 






(A.i) 



which are meant to approximate J^ f{x) dx for continuous functions / on some finite 
interval [a, b] c IR. We define the norm of a quadrature rule by 



\Q\ 



in 



Convergence of a sequence of quadrature rules is characterized by the following 
theorem of Polya (Davis and Rabinowitz 1984, p. 130). 

Theorem A.i. A sequence Q„ of quadrature rules converges for continuous functions, 



lim Qn{f) 



rb 

/ f{x)dx {feC[a,b]), 

J a 



^^A table can be obtained from the author upon request. Sasamoto (2005, Fig. 2) shows a plot (which 
differs by a scaling factor of two in both the function value and the time t) of the closely related function 

gi(t) = Y/var(yli(0-.4i(0))/2= Y^var(yli(0))-cov(yli(t),^i(0)) 
— without, however, commenting on either the numerical procedure used or on the accuracy obtained. 



ON THE NUMERICAL EVALUATION OF FREDHOLM DETERMINANTS 



37 



if and only if the sequence \\Qn\\ of norms is bounded by some stability constant A and if 



lim Qn{x') 



x^dx (fc = 0,1,2,... )■ 



(A.2) 



If the weights are all positive, then (A. 2) already implies the boundedness 0/ 1| Q„ [| = Q« (1). 

A quadrature rule Q is of order v ^ 1, if it is exact for all pol5moinials of degree 
at most V — 1. Using results from the theory polynomial best approximation one can 
prove quite strong error estimates (Davis and Rabinowitz 1984, §4.8). 

Theorem A. 2. If f E C*^^^'^ [a, b], then for each quadrature rule Q of order v ^ k with 
positive weights there holds the error estimate 

Q(/)- [''f{x)dx ^c,(l;-a)'^+iv-'^||/W||,»,„,„, 

with a constant^^ cj- depending only on k. 

If f is bounded analytic in the ellipse £p with foci at a, b and semiaxes of lengths s > cr 
such that 



S + (7 

s — a' 



then for each quadrature rule Q of order v with positive weights there holds the error estimate 



Qif)- / f{x)dx 



< 



A{b - a)p' 



Quadrature rules in two and more dimensions. For the n-dimensional integral 

f{ti,...,t„)dti ■ ■ ■ dt„ 

a,b]" 

we consider the product quadrature rule Q" that is induced by an one dimensional 
quadrature rule Q of the form (A.i), namely 






(A.3) 



We introduce some further notation for two classes of functions /. First, for / G 
C '^[[a,b]"), we define the seminorm 



! = 1 



(A.4) 



Second, if / G C{[a,b]"') is sectional analytic — that is, analytic with respect to each 
variable f, while the other variables are fixed in [a,b] — in the ellipse £„ (defined in 
Theorem A.2), and if / is uniformly bounded there, w^e call / to be of class Cp with 
norm 

n 

ll/llcp = E max \\f{tv---Ji-irJi+i,---Jn)\\L^{£„)- (A.5) 

/=1 {h,...,ti-i,ti+i,...,tn)e[a,b]"-^ 



^^Taking Jackson's inequality as given in Cheney (1998, p. 147), Cj^ = 7.[ne I ^)^ I \/lnk will do the job. 
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By a straightforward reduction argument (Davis and Rabinowitz 1984, p. 361) to 
the quadrature errors of the one-dimensional coordinate sections of /, Theorems A.i 
and A. 2 can now be generalized to n dimensions. 

Theorem A.3. If a sequence of quadrature rules converges for continuous functions, 
then the same holds for the induced n-dimensional product rules. 

If f E C^^^'^{[a,b]"), then for each one-dimensional quadrature rule Q of order v ^ k 
with positive weights there holds the error estimate 



Q"!/)-/ fih t„)dti---dt„ 

J[a,b]" 



^c,{b-ar+'v-'\f\,, 



with the same constant cj- depending only on k as in Theorem A. 2. 

Iff e C{[a,b]") is of class Cp, then for each one-dimensional quadrature rule Q of order 
V with positive weights there holds the error estimate 



Q"{f)- I fih t„)dtr---dt„ 

J\a,b]" 



Aib-aTp^ 



Notes on Gauss-Legendre and Curtis-Clenshaw quadrature. Arguably, the most in- 
teresting families of quadrature rules, with positive w^eights, are the Clenshaw-Curtis 
and Gauss-Legendre rules. With m points, the first is of order v = m, the second of or- 
der V = 2m. Thus, Theorems A.i and A. 2 apply The cost of computing the weights 
and points of Clenshaw-Curtis is 0{mlogm) using FFT, that of Gauss-Legendre 
is 0{m^) using the Golub-Welsh algorithm; for details see (Waldvogel 2006) and 
(Trefethen 2008). The latter paper studies in depth the reasons w^hy the Clenshaw- 
Curtis rule, despite having only half the order, performs essentially as w^ell as Gauss- 
Legendre for most integrands. To facilitate reproducibility we offer the Matlab code 
(which is just a minor variation of the code given in the papers mentioned above) 
that has been used in our numerical experiments: 

function [w,c] = ClenshawCurtis(a,b,m) 

m = m- 1 ; 

c = cos( (0:m) *pi/m) ; 

M = [l:2:m-l]'; 1 = length(M); n = m-1; 

vO = [2./M./(M-2); 1/M(end); zeros(n,l)]; 

v2 = -v0(l:end-l)-v0(end:-l:2); 

gO = -ones(m,l); gO(l+l)=gO(l+l)+m; gO(l+n)=gO(l+n)+m; 

g = g0/(m-2+mod(m,2)); 

w = ifft(v2+g); w(m+l) = w(l) ; 

c = ((l-c)/2*a+(l+c)/2*b)'; 

w = ((b-a)*w/2)'; 

for Clenshaw-Curtis; and 

function [w,c] = GaussLegendre(a,b,m) 

k = l:m-l; beta = k. /sqrt ( (2*k-l) . * (2*k+l) ) ; 

T = diag(beta,-l) + diag(beta, 1) ; 

[V,L] = eig(T); 

c = (diag(L)+l)/2; c = (l-c)+a+c+b; 

w = (b-a)*V(l,:) ."2; 
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for Gauss-Legendre, respectively. Note, however, that the code for Gauss-Legendre 
is, unfortunately, suboptimal in requiring O(m^) rather than O(m^) operations, since 
it establishes the full matrix V of eigenvectors of the Jacobi matrix T instead of 
directly calculating just their first components V(l, :) as in the fully fledged Golub- 
Welsh algorithm. Even then, there may well be more accurate, and more efficient, 
alternatives of computing the points and weights of Gauss-Legendre quadrature, see 
the discussions in Laurie (2001, §2) and Sw^arztrauber (2002, §4) and the literature 
cited therein. 

A.2. Determinantal bounds. In Section 6, for a continuous kernel X e C([fl, f^]^) 
of an integral operator, w^e need some bounds on the derivatives of the induced 
M-dimensional function 

K„(fi,...f„) = det(K(fp,f,))'^,^^^. 
To this end, if K e C'^~^'^([fl, V^) we define the norm 



|K||,, = max ||a'i3^2^||i,co. 



Lemma A.4. 1/K e C{[a,'bf), then K„ e C{[a,b]") with 



\\Kn\\,^^n"/'\\K 



n/2\\ jx||« 



(A.6) 



(A.7) 



IfK e C^ ^''^{[a,h]^), then Kn E C'^^^'^{[a,b]") with the seminorm (defined in (A.4)) 

|X„|fc^2'^M("+2)/2||K||». (A.8) 

IfK is bounded analytic on £p x £p (with the ellipse £p defined in Theorem A.2), then Kn is 
of class Cp (defined in (A.^)) and satisfies 



\\V W < »i(«+2)/2||fc'||" 



L~(£px£p)- 



(A.9) 



Proof. Using the multilinearity of the determinant we have 

K(ii,si) ■■■ K{ti,Sj) ■■■ K(fi,s„) 



3ff as} 



X(f,-,si) 



K{ti,Sj) 



K{t„,si) ■■■ K{t„,Sj) 



K{ti,s„) 
K{t„,s„j 



X(fi,si) 



a^K(fi,s^) 



K{ti,s„) 



d\K{ti,si) 



d\d[K{ti,Sj) 



d\K{ti,Sn) 



K{tn,Si) 



d'^K{tn,Sj) 



K{tn,S„) 
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which is, by Hadamard's inequality (Meyer 2000, p. 469)/4 bounded by the expres- 
sion (see also Lax (2002, p. 262)) 



m"''^ I max 



Now, with 



d'^K„{h,...,tn) = t. 



k\ d- 



k-l 3/ 



h\^Jdt]-'ds] 



Pld'^K\\i_ 

K(fi,si) 
K{t„,si) 



K{ti,s„) 
K{t„,s„) 



Si—ti,...,Su—tn 



we thus get 

\\:i^v II 



< 



n"'^ I max 



\d[d'^K\\ 



I'^n"'^ [ max p^d^^K\\i^^ 
yi+j^k 



1=0 \'/ \'+/<*^ 

This proves the asserted bounds (A. 7) and (A. 8) with k = and k ^ 1, respectively. 
The class Cp bound (A. 9) follows analogously to the case A: = 0. D 

A.3. Properties of a certain function used in Theorem 6.2. The power series 



~ „("+2)/2 

0(z)=E Z" 

defines an entire function on C (as the following lemma readily implies). 
Lemma A.5. Let Y be the entire function given by the expression 



(A.io) 



Y(z) = l + -^ze- 
If X > 0, then the series ^(x) is enclosed by:^^ 



^ -'/^ll + eri{-^ 



xT(xV2e) ^ cl)(x) ^ xY(xV2e). 



Proof. For x > we have 



cD(x) = X r ^x"-\ 

^ ' ^, Tin) 



1 r(«)' 

By Stirling's formula and monotonicity w^e get for n ^ 1 



,n/2 



< 



^1; 



71 r((n + l)/2)(V2e)"-i 
in fact, the upper bound is obtained for n — 1 and the low^er bound for n — > 00. Thus, 
by observing 

n=\ ^ ' 

we get the asserted enclosure. D 



^This inequality, discovered by Hadamard in 1893, ^^s already of fundamental importance to Fred- 
holm's original theory (Fredholm 1900, p. 41). 

^SNote the sharpness of this enclosure: \fejn = 0.93019 ■ ■ ■ . 
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